Descarga de datos de GBIF

Este ejercicio estará enfocado en descargar datos del repositorio GBIF y usar algunas herramientas sencillas para la limpieza de datos.

##Instalar los paquetes 

#install.packages("sf")
#install.packages("sp")
#install.packages("raster")
#install.packages("dismo")
#install.packages("spThin")
#install.packages("ecospat")
#install.packages("geodata")
#install.packages("sdm")
#install.packages("usdm")
#install.packages("biomod2")
#install.packages("rasterVis")




#Activar los paquetes 
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(sp)
library(raster)
library(dismo)
library(spThin)
## Loading required package: spam
## Spam version 2.10-0 (2023-10-23) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
## 
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
## 
##     backsolve, forwardsolve
## Loading required package: grid
## Loading required package: fields
## Loading required package: viridisLite
## 
## Try help(fields) to get started.
## Loading required package: knitr
library(ecospat)
library(geodata)
## Loading required package: terra
## terra 1.7.83
## 
## Attaching package: 'terra'
## The following object is masked from 'package:knitr':
## 
##     spin
## The following object is masked from 'package:fields':
## 
##     describe
## The following object is masked from 'package:grid':
## 
##     depth
## 
## Attaching package: 'geodata'
## The following object is masked from 'package:fields':
## 
##     world
library(sdm)
## sdm 1.2-46 (2024-07-17)
library(usdm)
library(rgdal)
## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: /usr/share/proj
## Linking to sp version:2.1-4
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## 
## Attaching package: 'rgdal'
## The following object is masked from 'package:terra':
## 
##     project
library(biomod2)
## biomod2 4.2-5-2 loaded.
##  /!\ New set up for modeling options. We apologize for the trouble ^[*.*]^
## Loading required package: nnet
## Loading required package: rpart
## Loading required package: mda
## Loading required package: class
## Loaded mda 0.5-4
## Loading required package: gam
## Loading required package: splines
## Loading required package: foreach
## Loaded gam 1.22-5
## Loading required package: mgcv
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:usdm':
## 
##     Variogram
## The following object is masked from 'package:raster':
## 
##     getData
## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.
## 
## Attaching package: 'mgcv'
## The following objects are masked from 'package:gam':
## 
##     gam, gam.control, gam.fit, s
## The following object is masked from 'package:nnet':
## 
##     multinom
## Loading required package: gbm
## Loaded gbm 2.2.2
## This version of gbm is no longer under development. Consider transitioning to gbm3, https://github.com/gbm-developers/gbm3
## Loading required package: earth
## Loading required package: Formula
## Loading required package: plotmo
## Loading required package: plotrix
## 
## Attaching package: 'plotrix'
## The following object is masked from 'package:terra':
## 
##     rescale
## The following object is masked from 'package:fields':
## 
##     color.scale
## Loading required package: maxnet
## Loading required package: randomForest
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
## Loading required package: xgboost
library(rasterVis)
## Loading required package: lattice
## 
## Attaching package: 'rasterVis'
## The following object is masked from 'package:gam':
## 
##     gplot
#Establacemos el directorio de trabajo 
setwd("/home2/taller_R/")
getwd()
## [1] "/home2/taller_R"

Buscaremos datos de presencia para una especie de lagartija del desierto (Phrynosoma modestum)

sp.gbif <- dismo::gbif(genus="Phrynosoma", species="modestum", geo=TRUE, removeZeros=TRUE, download=TRUE)
## 5012 records found
## 0-300-600-900-1200-1500-1800-2100-2400-2700-3000-3300-3600-3900-4200-4500-4800-5012 records downloaded

Vamos a indexar solos los datos de las coordenadas geográficas. Debemos saber en que posición están.

names(sp.gbif)
##   [1] "acceptedNameUsage"              "acceptedScientificName"        
##   [3] "acceptedTaxonKey"               "accessRights"                  
##   [5] "adm1"                           "adm2"                          
##   [7] "associatedOccurrences"          "associatedReferences"          
##   [9] "associatedSequences"            "associatedTaxa"                
##  [11] "basisOfRecord"                  "bibliographicCitation"         
##  [13] "catalogNumber"                  "class"                         
##  [15] "classKey"                       "cloc"                          
##  [17] "collectionCode"                 "collectionID"                  
##  [19] "collectionKey"                  "continent"                     
##  [21] "coordinatePrecision"            "coordinateUncertaintyInMeters" 
##  [23] "country"                        "crawlId"                       
##  [25] "datasetID"                      "datasetKey"                    
##  [27] "datasetName"                    "dateIdentified"                
##  [29] "day"                            "depth"                         
##  [31] "depthAccuracy"                  "disposition"                   
##  [33] "distanceFromCentroidInMeters"   "dynamicProperties"             
##  [35] "earliestAgeOrLowestStage"       "earliestEonOrLowestEonothem"   
##  [37] "earliestEpochOrLowestSeries"    "earliestPeriodOrLowestSystem"  
##  [39] "elevation"                      "elevationAccuracy"             
##  [41] "endDayOfYear"                   "establishmentMeans"            
##  [43] "eventDate"                      "eventRemarks"                  
##  [45] "eventTime"                      "family"                        
##  [47] "familyKey"                      "fieldNotes"                    
##  [49] "fieldNumber"                    "footprintSRS"                  
##  [51] "footprintWKT"                   "formation"                     
##  [53] "fullCountry"                    "gbifID"                        
##  [55] "gbifRegion"                     "genericName"                   
##  [57] "genus"                          "genusKey"                      
##  [59] "geodeticDatum"                  "georeferencedBy"               
##  [61] "georeferencedDate"              "georeferenceProtocol"          
##  [63] "georeferenceRemarks"            "georeferenceSources"           
##  [65] "georeferenceVerificationStatus" "habitat"                       
##  [67] "higherClassification"           "higherGeography"               
##  [69] "higherGeographyID"              "hostingOrganizationKey"        
##  [71] "http://unknown.org/captive"     "http://unknown.org/nick"       
##  [73] "http://unknown.org/recordID"    "identificationID"              
##  [75] "identificationQualifier"        "identificationRemarks"         
##  [77] "identifiedBy"                   "identifier"                    
##  [79] "individualCount"                "informationWithheld"           
##  [81] "installationKey"                "institutionCode"               
##  [83] "institutionID"                  "institutionKey"                
##  [85] "isInCluster"                    "ISO2"                          
##  [87] "isSequenced"                    "iucnRedListCategory"           
##  [89] "key"                            "kingdom"                       
##  [91] "kingdomKey"                     "language"                      
##  [93] "lastCrawled"                    "lastInterpreted"               
##  [95] "lastParsed"                     "lat"                           
##  [97] "latestAgeOrHighestStage"        "latestEpochOrHighestSeries"    
##  [99] "latestPeriodOrHighestSystem"    "license"                       
## [101] "lifeStage"                      "locality"                      
## [103] "locationAccordingTo"            "locationID"                    
## [105] "locationRemarks"                "lon"                           
## [107] "materialSampleID"               "modified"                      
## [109] "month"                          "municipality"                  
## [111] "nameAccordingTo"                "namePublishedIn"               
## [113] "namePublishedInYear"            "nomenclaturalCode"             
## [115] "occurrenceID"                   "occurrenceRemarks"             
## [117] "occurrenceStatus"               "organismID"                    
## [119] "organismName"                   "organismQuantity"              
## [121] "organismQuantityType"           "otherCatalogNumbers"           
## [123] "ownerInstitutionCode"           "parentNameUsage"               
## [125] "phylum"                         "phylumKey"                     
## [127] "preparations"                   "previousIdentifications"       
## [129] "projectId"                      "protocol"                      
## [131] "publishedByGbifRegion"          "publishingCountry"             
## [133] "publishingOrgKey"               "recordedBy"                    
## [135] "recordNumber"                   "references"                    
## [137] "rights"                         "rightsHolder"                  
## [139] "samplingProtocol"               "scientificName"                
## [141] "scientificNameID"               "sex"                           
## [143] "species"                        "speciesKey"                    
## [145] "specificEpithet"                "startDayOfYear"                
## [147] "taxonConceptID"                 "taxonID"                       
## [149] "taxonKey"                       "taxonomicStatus"               
## [151] "taxonRank"                      "taxonRemarks"                  
## [153] "type"                           "typeStatus"                    
## [155] "typifiedName"                   "verbatimCoordinateSystem"      
## [157] "verbatimElevation"              "verbatimEventDate"             
## [159] "verbatimLocality"               "verbatimSRS"                   
## [161] "vernacularName"                 "year"                          
## [163] "downloadDate"
sp.coords <- na.omit(sp.gbif[,c(106, 96, 140)]) # lon, lat, species
head(sp.coords)
##         lon      lat                   scientificName
## 1 -105.4433 32.23785 Phrynosoma modestum Girard, 1852
## 2 -105.4444 32.23727 Phrynosoma modestum Girard, 1852
## 3 -106.5001 32.86959 Phrynosoma modestum Girard, 1852
## 4 -107.4712 32.96034 Phrynosoma modestum Girard, 1852
## 5 -103.6514 29.57966 Phrynosoma modestum Girard, 1852
## 6 -106.7387 32.38824 Phrynosoma modestum Girard, 1852

Graficamos los puntos de las observaciones

plot(sp.coords$lon, sp.coords$lat) 

Aquí vemos que hay muchos puntos muy agrupados. Este tipo de datos puede influir negativamente en la calibración del modelo de nicho y la predicción geográfica. Vamos a tratar de reducir ese efecto adelgazando la muestra de puntos. Lo primero que hacemos es tratar de eliminar puntos a una distancia menor de 25 km. Como tenemos un tamaño de muestra grande (> 3000 registros) podemos hacer esto sin mucho problema. Usaremos el paquete spThin para esto.

sp.thin <-thin(sp.coords, lat.col="lat", long.col="lon", spec.col="scientificName", 
               thin.par=25, reps=1, locs.thinned.list.return=TRUE, write.files=TRUE,
               max.files=1, out.dir=".", out.base = "modestum.train.25km", write.log.file=TRUE, log.file="modestum.log.txt",
               verbose=TRUE)
## ********************************************** 
##  Beginning Spatial Thinning.
##  Script Started at: Mon Oct 28 15:18:34 2024
## Warning in thin(sp.coords, lat.col = "lat", long.col = "lon", spec.col =
## "scientificName", : There appear to be more than one species name in this *.csv
## file.
## Warning in thin(sp.coords, lat.col = "lat", long.col = "lon", spec.col =
## "scientificName", : Only using species name: Phrynosoma modestum Girard, 1852
## 
## lat.long.thin.count
## 455 
##   1 
## [1] "Maximum number of records after thinning: 455"
## [1] "Number of data.frames with max records: 1"
## [1] "Writing new *.csv files"
## Warning in thin(sp.coords, lat.col = "lat", long.col = "lon", spec.col =
## "scientificName", : Writing new *.csv files to output directory: .
## [1] "Writing file: ./modestum.train.25km_thin1.csv"

Ahora vamos primero a identificar que tipo de objeto se acaba de generar:

class(sp.thin)
## [1] "list"
str(sp.thin)
## List of 1
##  $ :'data.frame':    455 obs. of  2 variables:
##   ..$ Longitude: num [1:455] -107 -108 -104 -105 -109 ...
##   ..$ Latitude : num [1:455] 32.9 32.6 30.2 35.3 32.3 ...
head(sp.thin[[1]])
##    Longitude Latitude
## 3  -106.5001 32.86959
## 8  -107.9779 32.56886
## 26 -103.5483 30.15385
## 31 -105.2019 35.30268
## 33 -109.1064 32.31848
## 56 -108.8535 31.83653

Tenemos que tener en cuenta como vienen organizados los datos para luego manipularlos como queramos. Vamos a graficar los puntos adelgazados a una distancia no menor de 25 km.

plot(sp.thin[[1]]$Longitude, sp.thin[[1]]$Latitude)

Si comparamos con el objeto que contenia todos los puntos podemos ver que hemos reducido el sobre-muestreo en algunas regiones y esto nos facilitará la calibración y posterior predicción de nuestro modelo de nicho. Vamos a ver donde caen nuestros puntos dentro de un mapa.

Vamos a convertir nuestros puntos en un objeto espacial.

coords <- sp.thin[[1]]
coordinates(coords)<- ~Longitude+Latitude
plot(coords)

Si solo nos interesan los puntos que se distribuyen en México, podríamos excluir los puntos que caen por fuera del país. Vamos a hacerlo. Primero necesitamos cargar el mapa de México.

mx <- readOGR(dsn="./input/Poligono_Mexico", layer="gadm36_MEX_1")
## Warning in readOGR(dsn = "./input/Poligono_Mexico", layer = "gadm36_MEX_1"):
## OGR support is provided by the sf and terra packages among others
## Warning in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv =
## use_iconv, : OGR support is provided by the sf and terra packages among others
## Warning in ogrFIDs(dsn = dsn, layer = layer): OGR support is provided by the sf
## and terra packages among others
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : OGR support is provided by the sf and terra packages among others
## Warning in ogrListLayers(dsn): OGR support is provided by the sf and terra
## packages among others
## Warning in ogrFIDs(dsn = dsn, layer = layer): OGR support is provided by the sf
## and terra packages among others
## OGR data source with driver: ESRI Shapefile 
## Source: "/home2/taller_R/input/Poligono_Mexico", layer: "gadm36_MEX_1"
## with 32 features
## It has 10 fields
plot(mx)

Ahora le ponemos los puntos encima al mapa

plot(coords)
plot(mx, add=TRUE)

Con un clip vamos a quitar los puntos que caen por fuera. Pero primero definimos el sistema de coordenadas para nuestros puntos

mx
## class       : SpatialPolygonsDataFrame 
## features    : 32 
## extent      : -118.3689, -86.71014, 14.53292, 32.71804  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 10
## names       : GID_0, NAME_0,   GID_1,         NAME_1, VARNAME_1, NL_NAME_1,           TYPE_1,        ENGTYPE_1, CC_1, HASC_1 
## min values  :   MEX, Mexico, MEX.1_1, Aguascalientes,        NA,        NA, Distrito Federal, Federal District,   NA,  MX.AG 
## max values  :   MEX, Mexico, MEX.9_1,      Zacatecas,        NA,        NA,           Estado,            State,   NA,  MX.ZA
crdref <- CRS("+proj=longlat +datum=WGS84 +no_defs")
projection(coords)<-crdref 
coords
## class       : SpatialPoints 
## features    : 455 
## extent      : -118.1731, -95.86836, 21.20054, 44.71625  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs
coords.clip <- coords[mx, ]
coords.clip
## class       : SpatialPoints 
## features    : 195 
## extent      : -112.6083, -97.70008, 21.20054, 31.775  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs

Graficamos

plot(coords.clip)
plot(mx, add=TRUE)

Ahora tenemos nuestros puntos únicamente para México y los vamos a convertir en un objeto tipo data.frame con el fin de guardarlos en formato CSV y explorarlos con excel.

coords.clipDF <- as.data.frame(coords.clip)
head(coords.clipDF)
##     Longitude Latitude
## 119 -101.6513 23.88682
## 143 -101.2116 29.33872
## 187 -101.0510 22.76873
## 201 -102.7829 28.57397
## 309 -107.1414 31.54090
## 324 -102.0906 23.06960
write.csv(coords.clipDF, file="./output/Coordenadas_GBIF_Mexico.csv")
coordsMx <- read.csv("./output/Coordenadas_GBIF_Mexico.csv")
head(coordsMx)
##     X Longitude Latitude
## 1 119 -101.6513 23.88682
## 2 143 -101.2116 29.33872
## 3 187 -101.0510 22.76873
## 4 201 -102.7829 28.57397
## 5 309 -107.1414 31.54090
## 6 324 -102.0906 23.06960

Modelos de nicho ecológico o de distribución geográfica

En este parte del ejercicio vamos a realizar un modelo de nicho ecolĂłgico sencillo usando el algoritmo bioclim.

coordsMx <- read.csv("./output/Coordenadas_GBIF_Mexico.csv")
head(coordsMx)
##     X Longitude Latitude
## 1 119 -101.6513 23.88682
## 2 143 -101.2116 29.33872
## 3 187 -101.0510 22.76873
## 4 201 -102.7829 28.57397
## 5 309 -107.1414 31.54090
## 6 324 -102.0906 23.06960

Le quitamos la primera columna, no la necesitamos

coordsMx <- coordsMx[,c(2:3)]
head(coordsMx)
##   Longitude Latitude
## 1 -101.6513 23.88682
## 2 -101.2116 29.33872
## 3 -101.0510 22.76873
## 4 -102.7829 28.57397
## 5 -107.1414 31.54090
## 6 -102.0906 23.06960

Subimos las capas climáticas

clim <- worldclim_global(var="bio",res=2.5,path="./")

Gráficamos la primera

plot(clim[[1]])

clim<-stack(clim)

Vamos a hacer un buffer a nuestros puntos. Este buffer será nuestra área M del Diagrama BAM y es el área que la especie hipóteticamente ha podido “explorar” durante su historia evolutiva (i.e., dispersión, colonización, extinción):

coordinates(coordsMx) <- ~Longitude+Latitude
coordsMx <- sf::st_as_sf(coordsMx)
buffer.sp <- st_buffer(coordsMx, byid=FALSE, id=NULL, dist =3.0, nQuadsegs=5, 
                       endCapStyle="ROUND", joinStyle="ROUND", mitreLimit=1.0,dissolve = TRUE)
buffer.sp <- st_union(buffer.sp)

Graficamos

plot(buffer.sp)
plot(coordsMx, add=TRUE)
plot(mx,add=TRUE)

Vamos a cortar las capas climáticas usando este buffer

buffer.sp<-as_Spatial(buffer.sp)

class(buffer.sp)
## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"
clim.m <- crop(clim, buffer.sp)
plot(clim.m[[1]])

Ahora vamos a calibrar el modelo usando el algoritmo Bioclim

# primero convertimos las coordenadas a dataframe
coordsDF <- st_coordinates(coordsMx)
coordsDF <- as.data.frame(coordsDF)
class(coordsDF)
## [1] "data.frame"
# ajustamos un modelo de bioclim a los datos con las capas climáticas
bioc.modesta <- bioclim(clim.m, coordsDF)

Vamos a ver las gráficas de las variables de respuesta

response(bioc.modesta, var="wc2.1_2.5m_bio_1")

response(bioc.modesta, var="wc2.1_2.5m_bio_12")

response(bioc.modesta, var="wc2.1_2.5m_bio_4")

plot(bioc.modesta, a=1, b=2, p=0.95)

Ahora vamos a generar la predicción geográfica:

pred.bioc <- predict(clim.m, bioc.modesta)

plot(pred.bioc)

Vamos a generar una predicción binaria (0 ó 1; presencia vs. ausencia). Usaremos el criterio de umbral del “minimum presence training”

coords.bioc <- extract(pred.bioc, coords)
head(coords.bioc)
## [1] 0.02051282 0.02051282 0.01025641         NA 0.03589744 0.03076923
a <- min(coords.bioc)
a
## [1] NA

Usaremos ese valor para generar la predicción binaria. Por encima de 0.0073 será presencia y por debajo, ausencia.

bin.bioc <- reclassify(pred.bioc,c(-Inf,a,0, a,Inf,1))
bin.bioc
## class      : RasterLayer 
## dimensions : 398, 502, 199796  (nrow, ncol, ncell)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -115.625, -94.70833, 18.20833, 34.79167  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 0.6410256  (min, max)

Graficamos

plot(bin.bioc)

Vamos a usar los registros completos de la distribuciĂłn de la especie para ver si tenemos un problema de truncaciĂłn de nicho

Hacemos otro buffer pero con todos los puntos adelgazados. Checamos la anchura para que todo quede conectado

coords<-as.data.frame(coords)
head(coords)
##    Longitude Latitude
## 3  -106.5001 32.86959
## 8  -107.9779 32.56886
## 26 -103.5483 30.15385
## 31 -105.2019 35.30268
## 33 -109.1064 32.31848
## 56 -108.8535 31.83653
write.csv(coords, file="./output/Coordenadas_GBIF.csv")
coords <- read.csv("./output/Coordenadas_GBIF.csv")
head(coords)
##    X Longitude Latitude
## 1  3 -106.5001 32.86959
## 2  8 -107.9779 32.56886
## 3 26 -103.5483 30.15385
## 4 31 -105.2019 35.30268
## 5 33 -109.1064 32.31848
## 6 56 -108.8535 31.83653
coords <- coords[,c(2:3)]
head(coords)
##   Longitude Latitude
## 1 -106.5001 32.86959
## 2 -107.9779 32.56886
## 3 -103.5483 30.15385
## 4 -105.2019 35.30268
## 5 -109.1064 32.31848
## 6 -108.8535 31.83653
coordinates(coords) <- ~Longitude+Latitude
head(coords)
## class       : SpatialPoints 
## features    : 1 
## extent      : -106.5001, -106.5001, 32.86959, 32.86959  (xmin, xmax, ymin, ymax)
## crs         : NA
coords <- sf::st_as_sf(coords)
buffer.sp2 <- st_buffer(coords, byid=FALSE, id=NULL, dist =6.0, nQuadsegs=5, 
                       endCapStyle="ROUND", joinStyle="ROUND", mitreLimit=1.0,dissolve = TRUE)
buffer.sp2 <- st_union(buffer.sp2)

Graficamos

plot(buffer.sp2)
plot(coords, add=TRUE)
plot(mx, add=TRUE)

buffer.sp2<-as_Spatial(buffer.sp2)

class(buffer.sp2)
## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"
clim.m2 <- crop(clim, buffer.sp2)
plot(clim.m2[[1]])

coordsDF2 <- st_coordinates(coords)
coordsDF2 <- as.data.frame(coordsDF2)
class(coordsDF2)
## [1] "data.frame"
bioc.modesta.v2 <- bioclim(clim.m2, coordsDF2)
pred.bioc.v2 <- predict(clim.m2, bioc.modesta.v2)
plot(pred.bioc.v2)

Comparemos las dos predicciones:

par(mfrow=c(1,2))
plot(pred.bioc, main="M puntos en México")
plot(pred.bioc.v2, main="M todos los puntos")

Ajustar las extensiones para ambos rasters

pred.bioc <- projectRaster(from=pred.bioc, to=pred.bioc.v2, crs=crs(pred.bioc.v2), method="bilinear")

par(mfrow=c(1,2))
plot(pred.bioc, main="M puntos en México")
plot(pred.bioc.v2, main="M todos los puntos")

Aquí hemos calibrado un modelo de nicho ecológico usando todos los datos disponibles y por lo tanto no tenemos forma de evaluar que tanto error puede tener nuestra predicción geográfica. Vamos a partir los datos en dos conjuntos. Uno de calibración y uno de validación:

La partición de los datos será aleatoria. Usaremos la función kfold para esto. Usaremos 5 grupos

group <- kfold(coordsDF2, 5)
group
##   [1] 3 1 1 3 4 2 5 1 5 2 1 3 5 4 5 2 4 2 1 4 1 4 4 2 2 2 3 4 4 1 4 2 1 3 5 5 3
##  [38] 2 3 1 3 4 3 3 1 3 4 4 4 3 1 2 4 3 2 4 2 5 4 3 3 1 3 4 4 5 1 4 1 3 5 5 5 5
##  [75] 3 3 3 5 5 5 2 2 2 4 4 5 1 1 3 4 2 3 3 1 2 3 5 2 5 3 5 5 2 5 2 1 1 5 4 3 3
## [112] 5 5 3 1 1 2 5 1 4 4 2 2 5 2 1 1 1 1 5 3 4 2 2 2 4 4 4 2 2 4 3 3 1 3 4 2 2
## [149] 2 2 3 2 1 3 2 5 1 4 1 2 1 4 1 3 2 4 5 4 5 3 4 2 4 5 3 1 3 5 3 5 1 5 2 5 4
## [186] 5 4 2 4 4 3 4 5 4 5 3 1 3 1 3 4 1 1 3 2 5 5 3 3 4 4 1 5 2 2 4 2 3 5 3 3 3
## [223] 1 1 3 5 5 4 2 5 2 1 3 2 4 5 2 3 1 1 1 4 2 2 3 2 4 3 2 1 4 3 4 3 2 2 3 4 3
## [260] 5 1 3 1 1 1 4 2 2 1 1 4 4 4 4 5 4 1 4 1 5 4 1 1 4 4 4 3 5 2 5 2 4 2 3 5 2
## [297] 2 3 5 1 4 3 5 1 4 2 5 2 3 4 2 5 1 5 4 2 3 2 3 1 1 5 5 5 3 4 2 1 5 1 5 1 2
## [334] 2 3 2 2 5 5 1 5 1 2 5 3 3 3 2 5 2 5 1 1 1 5 1 2 1 1 3 4 5 2 3 1 2 3 5 5 2
## [371] 2 3 4 5 1 5 4 2 3 1 4 4 2 1 1 4 5 4 3 5 3 5 1 5 3 4 3 4 5 3 3 1 2 5 4 2 5
## [408] 1 1 4 3 4 3 4 1 5 5 5 1 2 3 2 4 2 3 3 2 2 5 4 1 1 2 1 5 1 4 4 1 2 3 5 2 5
## [445] 4 3 3 5 4 1 5 1 1 5 4
length(group)
## [1] 455
pres_train <- coordsDF[group != 1, ]
pres_test <- coordsDF[group == 1, ]

head(pres_train)
##           X        Y
## 1 -101.6513 23.88682
## 4 -102.7829 28.57397
## 5 -107.1414 31.54090
## 6 -102.0906 23.06960
## 7 -101.5608 24.71512
## 9 -101.1581 25.28576
head(pres_test)
##            X        Y
## 2  -101.2116 29.33872
## 3  -101.0510 22.76873
## 8  -100.8064 22.63093
## 11 -106.2099 27.87126
## 19 -100.5416 26.22126
## 21 -100.7210 25.91715

Ya tenemos las presencias particionadas en dos conjuntos (calibración o “training” y validación o “testing”). Para validar nuestros modelos necesitamos las “ausencias”. No tenemos ausencias, así que vamos a generar unas ausencias aleatorias a partir de la geografía: a este tipo de ausencias se les conoce como “pseudo-ausencias”

ext <- extent(clim.m2)
backg <- randomPoints(clim.m2, n=1000, ext=ext)
colnames(backg) = c('Longitude', 'Latitude')
group <- kfold(backg, 5)
backg_train <- backg[group != 1, ]
backg_test <- backg[group == 1, ]

Vamos a graficar los puntos

plot(clim.m2[[1]])
points(backg_train, pch='-', cex=1, col='yellow')
points(backg_test, pch='-',  cex=1, col='black')
points(pres_train, pch= '+', col='green')
points(pres_test, pch='+', col='blue')

Ahora vamos a calibrar con el primer conjunto que extrajimos, pres_train:

bc <- bioclim(clim.m2, pres_train)
plot(bc, a=1, b=2, p=0.85)

Vamos a validar nuestro modelo de Bioclim

# Usaremos la funciĂłn evaluate de dismo.
eval.modesta <- evaluate(pres_test, backg_test, bc, clim.m2)
eval.modesta
## class          : ModelEvaluation 
## n presences    : 34 
## n absences     : 200 
## AUC            : 0.8564706 
## cor            : 0.5310529 
## max TPR+TNR at : 0.01853354

Vamos a calcular un umbral

tr <- threshold(eval.modesta, 'spec_sens')
tr
## [1] 0.01853354

Vamos a hacer la predicciĂłn en la geografĂ­a:

pred.bioc.v3 <- predict(clim.m2, bc, ext=ext, progress='')
pred.bioc.v3
## class      : RasterLayer 
## dimensions : 852, 823, 701196  (nrow, ncol, ncell)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -124.1667, -89.875, 15.20833, 50.70833  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 0.621118  (min, max)

Graficamos los dos mapas de predicciones, el continuo y el discreto usando el valor de umbral que calculamos antes.

par(mfrow=c(1,2))
plot(pred.bioc.v3, main='Bioclim, valores crudos')
plot(pred.bioc.v3 > tr, main='presencia/ausencia')
points(pres_train, pch='+')

Ahora vamos a generar modelos con varios algoritmos y crear un ensamble. Usaremos un conjunto de datos de un ratón (Neotoma cinerea “bushy-tailed woodrat”) de la familia Cricetidae que se distribuye en Estados Unidos y Cánada. Usaremos el paquete de R “sdm” (Naimi & Aráujo 2016; https://github.com/babaknaimi/sdm) que es bastante versátil.

datos <- read.csv("Neotoma.csv")
plot(datos$decimalLongitude,datos$decimalLatitude)

bios <- stack(list.files(path="./ClimaActual/", pattern = "*.tif", full.names = T))
plot(bios[[1]])

pts <- datos[,c(2,3)]
selected <- sample(1:nrow(pts), nrow(pts) * 0.5)
train <- pts[selected,] # this is the selection to be used for model training
test <- pts[-selected,] # this is the opposite of the selection which will be used for model testing
pts.all <- rbind(train,test)
# convertir el primer raster de los bios en puntos 
bios.pts<-rasterToPoints(bios[[1]],fun=NULL, spatial=TRUE)
bios.all<-raster::extract(bios[[1]],bios.pts)
bios.all<-data.frame(coordinates(bios.pts),bios.all)

# generar pseudo-ausencias para los datos de calibración (50,000 o más)
presences_train<- train
dim(presences_train)
## [1] 216   2
pseudo_abs_train<- ecospat.rand.pseudoabsences(nbabsences=1000,
                                               glob=bios.all,colxyglob=1:2, colvar=1:2, presence=presences_train,
                                               colxypresence=1:2, mindist=0.1)

pseudo_abs_train<- pseudo_abs_train[,1:2]
presences_train["Rep1"] <- 1
presences_train <- plyr::rename(presences_train, replace=c("decimalLongitude"="x", "decimalLatitude"="y"))
pseudo_abs_train["Rep1"] <- 0

# Combino las presencias y las pseudo-absences del conjunto de calibraciĂłn
DataSpeciesTrain.spp <-rbind(presences_train, pseudo_abs_train)

# generar pseudo-ausencias para los datos de validaciĂłn
presences_test <- test
presences_test["Rep1"] <- 1
presences_test <- plyr::rename(presences_test, replace=c("decimalLongitude"="x", "decimalLatitude"="y"))
dim(presences_test)
## [1] 217   3
pseudo_abs_test<- ecospat.rand.pseudoabsences(nbabsences=1000,
                                              glob=bios.all,colxyglob=1:2, colvar=1:2, presence=presences_test,
                                              colxypresence=1:2, mindist=0.1)

pseudo_abs_test<-pseudo_abs_test[,1:2]
pseudo_abs_test["Rep1"] <- 0

# Combino las presencias y las pseudo-absences del conjunto de validaciĂłn
DataSpeciesTest.spp <- rbind(presences_test, pseudo_abs_test)


# Convierto los dataframes en objectos espaciales
DataSpecies <- SpatialPointsDataFrame(coords =DataSpeciesTrain.spp[,c(1,2)], data = DataSpeciesTrain.spp,
                                      proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
TestSpecies <- SpatialPointsDataFrame(coords =DataSpeciesTest.spp[,c(1,2)], data = DataSpeciesTest.spp,
                                      proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
DataSpecies <- DataSpecies[,c(-1,-2)]
TestSpecies <- TestSpecies[,c(-1,-2)]

Voy a examinar la colinealidad entre las variables usando el VIF (https://en.wikipedia.org/wiki/Variance_inflation_factor)

biosDF <- as.data.frame(bios)
# La funciĂłn vifcor identifica las variables que deben ser excluidas
v1 <- usdm::vifcor(biosDF, th=0.7)
bios.red <- usdm::exclude(bios, v1)
bios.red
## class      : RasterStack 
## dimensions : 1363, 2741, 3735983, 3  (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667  (x, y)
## extent     : -166.8333, -52.625, 14.54167, 71.33333  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      : bio_05, bio_12, bio_15 
## min values :    -27,     35,      5 
## max values :    456,   5152,    138

Vamos a ajustar los datos a los requirimientos especĂ­ficos del paquete sdm usando todas las bios

d1 <- sdmData(formula=Rep1~., train=DataSpecies, test=TestSpecies, predictors=bios)

d1
## class                                 : sdmdata 
## =========================================================== 
## number of species                     :  1 
## species names                         :  Rep1 
## number of features                    :  5 
## feature names                         :  bio_01, bio_05, bio_06, ... 
## type                                  :  Presence-Absence 
## has independet test data?             :  TRUE 
## number of records                     :  train-> 1214; test-> 1216 
## has Coordinates?                      :  TRUE

Estrategia de submuestreo: 70%-30% (calibración-validación) con 2 (haga 10 o más en su casa) replicas de bootstrapping. Vamos a probar maxent con default settings (ya vimos en clase qué consecuencias trae eso), random forest (los mejores en mi opinión), glm y bioclim. Se pueden ajustar más algoritmos. Vaya a la ayuda de la función tecleando ?sdm para ver más opciones.

m1 <- sdm(Rep1~., data=d1, methods=c('rf','glm','bioclim'),
          replicatin='sub', test.percent=30, n=2)

m1
## class                                 : sdmModels 
## ======================================================== 
## number of species                     :  1 
## number of modelling methods           :  3 
## names of modelling methods            :  rf, glm, bioclim 
## replicate.methods (data partitioning) :  subsampling 
## number of replicates (each method)    :  2 
## toral number of replicates per model  :  2 (per species) 
## test percentage (in subsampling)      :  30 
## ------------------------------------------
## model run success percentage (per species)  :
## ------------------------------------------
## method          Rep1             
## ---------------------- 
## rf         :        100   %
## glm        :        100   %
## bioclim    :        100   %
## 
## ###################################################################
## model Mean performance (per species), using independent test dataset:
## -------------------------------------------------------------------------------
## 
##  ## species   :  Rep1 
## =========================
## 
## methods    :     AUC     |     COR     |     TSS     |     Deviance 
## -------------------------------------------------------------------------
## rf         :     0.92    |     0.65    |     0.75    |     0.53     
## glm        :     0.88    |     0.56    |     0.62    |     0.64     
## bioclim    :     0.81    |     0.52    |     0.6     |     0.96

Vamos a ver las curvas ROC que evaluan el balance entre la especificidad y sensibilidad. Recuerden la clase de la curva ROC y matriz de confusiĂłn.

roc(m1)

Generar un ensamble con los mejores modelos usando el método de promedio ponderado (weighted averaging) con base en la métrica (TSS: True Statistics Skill)

e1 <- ensemble(m1, newdata=bios, filename='e1.img', setting=list(method='weighted',stat='TSS'))
## Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Message
## 1: NaN converted to INT_MAX.

Graficamos el modelo

plot(e1)
## Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Message
## 1: NaN converted to INT_MAX.
## Warning in x@ptr$readStart(): GDAL Message 1: NaN converted to INT_MAX.

Generamos un mapa binario usando un criterio de umbral de mínima área predicha y donde le asignamos un error de 10% a los datos de presencia que no usamos para calibrar.

xy <- train
coordinates(xy) <- ~decimalLongitude+decimalLatitude
mpa.e1 <- ecospat.mpa(e1, xy, perc = .9)
## Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Message
## 1: NaN converted to INT_MAX.
## Warning in x@ptr$extractCell(i - 1): GDAL Message 1: NaN converted to INT_MAX.
mpa.e1
##   10% 
## 0.301
##Cambiar e1 de rasterlayer a SpatRaster
e1<-rast(e1)
## Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Message
## 1: NaN converted to INT_MAX.
bin.e1 <- ecospat.binary.model(e1, mpa.e1)
## Warning in x@ptr$classify(as.vector(rcl), NCOL(rcl), right, include.lowest, :
## GDAL Message 1: NaN converted to INT_MAX.
plot(bin.e1)

Vamos a generar unas proyecciones a futuro.

# Subir bioclimáticos de escenarios futuros (RCP8.5 2070)
fut.bios <- stack(list.files(path="./CMIP5", pattern = "tif", full.names = T))

fut.bios <- projectRaster(from=fut.bios, to=bios[[1]], crs=fut.bios, method="bilinear")

# recortar las bios a la extensiĂłn
fut.bios.m <- raster::crop(fut.bios, bios[[1]])


# generar la proyecciĂłn futura sobre el ensamble
e1.fut <- ensemble(m1,newdata=fut.bios.m,filename='e1_fut.img',overwrite = T, setting=list(method='weighted',stat='TSS'))
## Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Message
## 1: NaN converted to INT_MAX.

Graficamos

plot(e1.fut)
## Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Message
## 1: NaN converted to INT_MAX.
## Warning in x@ptr$readStart(): GDAL Message 1: NaN converted to INT_MAX.

Generar el binario para el escenario futuro

e1.fut<-rast(e1.fut)
## Warning in new_CppObject_xp(fields$.module, fields$.pointer, ...): GDAL Message
## 1: NaN converted to INT_MAX.
bin.e1.fut <- ecospat.binary.model(e1.fut, mpa.e1)
## Warning in x@ptr$classify(as.vector(rcl), NCOL(rcl), right, include.lowest, :
## GDAL Message 1: NaN converted to INT_MAX.

Graficamos el mapa binario proyectado a futuro

plot(bin.e1.fut)

Ahora vamos a identificar áreas de estabilidad, ganancia y pérdida de áreas (i.e., número de pixeles). Ya tenemos los mapas binarios para el periodo actual y el futuro (0, 1). Usaremos una función del paquete Biomod2 (BIOMOD_RangeSize).

bin.e1<-raster(bin.e1)
bin.e1.fut<-raster(bin.e1.fut)
bin.e1 <- reclassify(bin.e1, c(-Inf, 0, 0, 0, Inf, 1))
bin.e1.fut <- reclassify(bin.e1.fut, c(-Inf, 0, 0, 0, Inf, 1))

cambio <- BIOMOD_RangeSize(bin.e1, bin.e1.fut)
## 
## -=-=-=-=-=-=-=-=-=-=-=-=-= Do Range Size Computation -=-=-=-=-=-=-=-=-=-=-=-=-=
## 
##  Comparing 'proj.current' and 'proj.future'. 
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
cambio
## $Compt.By.Models
##                    Loss Stable0 Stable1   Gain PercLoss PercGain
## ensemble_weighted 30411 1301587  130205 162337   18.934  101.071
##                   SpeciesRangeChange CurrentRangeSize FutureRangeSize.NoDisp
## ensemble_weighted             82.138           160616                 130205
##                   FutureRangeSize.FullDisp
## ensemble_weighted                   292542
## 
## $Diff.By.Pixel
## class       : SpatRaster 
## dimensions  : 1363, 2741, 1  (nrow, ncol, nlyr)
## resolution  : 0.04166667, 0.04166667  (x, y)
## extent      : -166.8333, -52.625, 14.54167, 71.33333  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
## source(s)   : memory
## name        : ensemble_weighted 
## min value   :                -2 
## max value   :                 1

Vamos a mostrar el mapa.

Valores de -2 son los pixeles que se proyectan por perderse para la especie en el futuro. Valores de -1 son los pixeles que se proyectan por estar estables (no cambian). Valores de 0 son los pixeles que no ocupados actualmente pero tampoco en el futuro. Valores de 1 son los pixeles no ocupados actualmente pero proyectados por estarlo en el futuro.

plot(cambio$Diff.By.Pixel)

Hay varios paquetes para visualizar rasters.

https://oscarperpinan.github.io/rastervis/

levelplot(cambio$Diff.By.Pixel, par.settings=GrTheme())

TAREA:

Como tarea les queda hacer las proyecciones futuras con cada modelo invididual y evaluar qué tanta incertidumbre hay en cada algoritmo.

Eso es todo por ahora!